require(geoR)
## Loading required package: geoR
## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------
data(s100)
summary(s100)
## Number of data points: 100
##
## Coordinates summary
## Coord.X Coord.Y
## min 0.005638006 0.01091027
## max 0.983920544 0.99124979
##
## Distance summary
## min max
## 0.007640962 1.278175109
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.1676955 0.2729882 1.1045936 0.9307179 1.6101707 2.8678969
##
## Other elements in the geodata object
## [1] "cov.model" "nugget" "cov.pars" "kappa" "lambda"
names(s100)
## $coords
## NULL
##
## $data
## [1] "data"
##
## $other
## [1] "cov.model" "nugget" "cov.pars" "kappa" "lambda"
var(s100$data)
## [1] 0.7322958
boxplot(s100$data)

## points.geodata()
points(s100)

points(s100, pt.div="equal")

points(s100, cex.min=.6, cex.max=.6)

points(s100, cex.min=.3, cex.max=3)

points(s100, cex.min=1, cex.max=1, pt.div="quart")

points(s100, pt.div="quart")

points(s100, pt.div=4)

points(s100, pt.div="dec")

points(s100, cex.min=1, cex.max=1, pt.div="quint")

points(s100, cex.min=1, cex.max=1, col="gray")

points(s100, cex.min=1, cex.max=1, col=rainbow(100))

points(s100, cex.min=1, cex.max=1, col=terrain.colors(100))

points(s100, pt.div="equal", pch.seq="+")

points(s100, data=sample(s100$data))

points(s100, data=exp(s100$data))

class(s100)
## [1] "geodata"
args(points.geodata)
## function (x, coords = x$coords, data = x$data, data.col = 1,
## borders = x$borders, pt.divide = c("data.proportional", "rank.proportional",
## "quintiles", "quartiles", "deciles", "equal"), lambda = 1,
## trend = "cte", abs.residuals = FALSE, weights.divide = "units.m",
## cex.min, cex.max, cex.var, pch.seq, col.seq, add.to.plot = FALSE,
## x.leg, y.leg = NULL, dig.leg = 2, round.quantiles = FALSE,
## permute = FALSE, ...)
## NULL
help(points.geodata)
data(parana)
summary(parana)
## Number of data points: 143
##
## Coordinates summary
## east north
## min 150.1220 70.3600
## max 768.5087 461.9681
##
## Distance summary
## min max
## 1.0000 619.4925
##
## Borders summary
## east north
## min 137.9873 46.7695
## max 798.6256 507.9295
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 162.7700 234.1900 269.9200 274.4106 318.2300 413.7000
##
## Other elements in the geodata object
## [1] "loci.paper"
borders <- parana$borders
points(parana, border=borders)

points(parana, border=borders, trend="1st")

points(parana, border=borders, trend="1st", pt.div="quart")

## plot.geodata()
plot(s100)

plot(parana)

plot(parana, bor=borders, low=T)

data(ca20)
summary(ca20)
## Number of data points: 178
##
## Coordinates summary
## east north
## min 4957 4829
## max 5961 5720
##
## Distance summary
## min max
## 43.01163 1138.11774
##
## Borders summary
## east north
## min 4920 4800
## max 5990 5800
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 21.00000 43.00000 50.50000 50.67978 58.00000 78.00000
##
## Covariates summary
## altitude area
## Min. :3.300 1: 13
## 1st Qu.:5.200 2: 49
## Median :5.650 3:116
## Mean :5.524
## 3rd Qu.:6.000
## Max. :6.600
##
## Other elements in the geodata object
## [1] "reg1" "reg2" "reg3"
names(ca20)
## $coords
## [1] "east" "north"
##
## $data
## [1] "data"
##
## $covariate
## [1] "altitude" "area"
##
## $borders
## [1] "borders"
##
## $other
## [1] "reg1" "reg2" "reg3"
plot(ca20)

points(ca20, cex.min=0.5, cex.max=0.5)
polygon(ca20$reg1)
polygon(ca20$reg2)
polygon(ca20$reg3)

plot(ca20)

plot(ca20,trend=~area)

plot(ca20,trend=~area+altitude)

ksat <- read.geodata("http://www.leg.ufpr.br/geoR/tutorials/datasets/Cruciani.dat", head=T, row.names=1)
summary(ksat)
## Number of data points: 32
##
## Coordinates summary
## x y
## min 2.5 0.8
## max 21.0 9.1
##
## Distance summary
## min max
## 0.509902 18.552897
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0200000 0.0922500 0.3050000 0.8720625 1.1425000 4.3000000
##
## Other elements in the geodata object
## [1] "call"
plot(ksat)

require(MASS)
## Loading required package: MASS
boxcox(ksat)

plot(ksat, lambda=0)

kbor <- read.table("http://www.leg.ufpr.br/geoR/tutorials/datasets/Cruciani.border", head=T, row.names=1)
kbor
## x y
## 1 0.0 0.0
## 2 0.7 1.6
## 3 2.2 3.3
## 4 5.0 5.7
## 5 8.7 7.8
## 6 13.9 10.4
## 7 19.8 13.0
## 8 21.0 10.5
## 9 22.4 5.3
## 10 22.5 0.0
plot(ksat, lambda=0, bor=kbor)

points(ksat, lambda=0, bor=kbor)

## Variog
s100.v <- variog(s100, max.d=1)
## variog: computing omnidirectional variogram
plot(s100.v)

parana.v <- variog(parana, max.dist=400, trend="1st")
## variog: computing omnidirectional variogram
plot(parana.v)

ca20.v <- variog(ca20, trend=~area, max.d=800)
## variog: computing omnidirectional variogram
plot(ca20.v)

ksat.v <- variog(ksat, max.d=10, lam=0)
## variog: computing omnidirectional variogram
plot(ksat.v)
